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m | ABSTRACT 



In this paper we numerically simulate some of the most critical physical processes in 
galaxy formation: The supernova feedback loop, in conjunction with gas dynamic pro- 
' cesses and gravitational condensations, plays a crucial role in determining how the 

OO . observable properties of galaxies arise within the context of a model for large-scale 

structure. Our treatment incorporates a multi-phase model of the interstellar medium 
and includes the effects of cooling, heating and metal enrichment by supernovae, and 
evaporation of cold clouds. The star formation happens inside the clouds of cold gas, 
, which are produced via thermal instability. In this paper we simulate the galaxy for- 

mation in standard biased Cold Dark Matter (CDM) model for a variety of parameters 
and for several resolutions. 

In our picture, supernova feedback regulates the evolution of the gas components and 
star formation. The efficiency of cold cloud evaporation by supernova strongly influences 
star formation rates. This feedback results in a steady rate of star formation in "large" 
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, galaxies (mass larger than (2 — 3) x 10 u Mq within 100 kpc radius) at a level of 

(1 — 10)Mq per year for z < 3 (Hq — 50 Km s _1 Mpc -1 ). Supernova feedback has an 
even stronger effect on the evolution of "dwarf" galaxies. Most of the dwarf galaxies in 
' our models have a small fraction of stars and extremely low luminosities: Mr > — 15 

for parent dark-halo masses M to t < (2 — 3) x 1O 1O M within a 50 kpc radius. 
The observational properties (colors, luminosities) of galaxies identified in the simula- 
tions are computed using a stellar population synthesis model. In the case of both large 
and small galaxies, the distribution of luminous matter (stars) is strongly biased with 
respect to the dark matter. For a range of parameter values and resolutions we find an 
approximate biasing measure of the form p\ um = (/Odm/133) 1 ' 7 , for overdensities exceed- 
ing about 1000. Deviations from this relation depend strongly on the environment. For 
halo masses exceeding 2 x 10 10 M Q , the dependence of the absolute visual magnitude 
My on the total mass can be approximated as My — —18.5 — 41og(Mtot/10 11 Mo ), 
with a scatter of less than 1/2 magnitude. 

Key words: Galaxies: evolution - Galaxies: formation - Hydrodynamics - Methods: 
numerical 



1 INTRODUCTION 

In the past, an important reason for the attractiveness of dark-matter models for large-scale structure has been their complete 
neglect of nongravitational forces, permitting rather far reaching conclusions to be drawn using relatively simple physical 
considerations combined with quite sophisticated mathematical and statistical methods. Improvements in computational 
resources now permit numerical simulations of large-scale structure to approach a resolution at which the neglect of non- 
gravitational effects, rather than the limited numerical resolution, is the principal source of uncertainties. 

In this paper, we study the formation and evolution of galaxies in the context of large-scale structure. We focus our 
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attention on the role of supernova feedback and hydrodynamics in combination with gravitational dynamics in determining 
the evolution of some of the most important observational properties of galaxies, including luminosities and colors, and the 
extent of their dependence on mass and environment. 

The formation and evolution of galaxies and the dynamics of the intergalactic medium are complicated processes which 
cannot be modelled purely from first principles and are not yet completely understood. One of the main challenges to 
modeling is the fact that a variety of effects over an exceptionally wide range of scales all come into play: Cosmological 
fluctuations on scales larger than tens of megaparsecs are important because these scales define the large-scale structure 
of the galaxy distribution - positions and sizes of clusters, superclusters, filaments, and voids. Many parameters of galaxies 
(like the morphological type, rate of collisions and merging) are known to depend on the position of a galaxy within the 
structures. Conditions for galaxy formation are different for galaxies inside groups and in the field, which brings in a scale of 
one megaparsec. Although the dark matter still dominates gravitationally on scales above 100 kpc, the baryonic component 
becomes more and more important as we descend from the megaparsec scale. The dissipation of energy by the hot gas definitely 
is an important factor on scales below 100 kpc. For the formation of the luminous component of galaxies, scales below 10 kpc 
are essential. These scales define the global dynamics of the galaxy - density and velocity profiles, spiral arms, and so on. 
In turn, as discussed in more detail below, the global galactic dynamics regulates the formation of cold gas clouds (in our 
galaxy - molecular clouds) inside of which stars are produced. The characteristic scale of the clouds is 10-100 pc. The clouds 
themselves have a very complicated internal structure. 

Thermal instability plays an important role in these processes: It has long been appreciated (Field 1965, Balbus 1986, 
Nulsen 1986 ) that gas subject to cooling (and perhaps heating) processes may become thermally unstable. In particular, 
for primordial abundances, pressure equilibrium, and in the absence of photoionization, it is known that thermal instabilities 
are possible for gas within a wide range of temperatures that are quite typical in a cosmological setting. For example, 
thermal instabilities are important in the theory of globular clusters (Fall & Rees, 1985). The situation is more complicated 
if photoionization is properly taken into account (Miicket & Kates, 1996), but the general picture remains that at densities 
and temperatures typical for a galactic halo, and with reasonable estimates for the ionizing flux, instabilities will evolve on a 
timescale that is sufficiently short compared to appropriate dynamic timescales such as a free fall time. 

One important consequence of thermal instability is that it causes the formation of cool clouds, whose temperature 
continues to drop until the cool gas becomes neutral and essentially stops losing energy. This occurs at a temperature of 
roughly 10 4 Jf. Now, if the heat gain due to conductivity is greater than the loss due to radiation, clouds will not persist. 
However, beyond a critical size, clouds tend to grow because conduction is not sufficient to evaporate them (Begelman & 
McKee, 1990). When the size of a cloud approaches the Jeans limit, it starts to form stars with all the ensuing complications. 
Eventually, the cloud presumably stops growing as a result of reheating due to stellar winds, supernova explosions, and 
radiation from young stars. Subsequent cooling of the hot gas could form new clouds. 

Cold clouds are almost certainly an essential ingredient in the dynamics of galaxies, particularly in star formation and 
the interstellar medium. Their importance has long been appreciated, and there have been many refinements in the multi- 
component picture of the interstellar medium (Oort, 1954; Spitzer, 1965; Field & Saslaw, 1965; Field, Goldsmith & Habing, 
1969; McKee & Cowie, 1977) Thermal instabilities and cloud formation are important in the theory of globular clusters (Fall 
& Rees, 1985). 

Another important component of galaxy formation and evolution is supernovae and young stars. Supernovae eject enor- 
mous amounts of energy into the surrounding gas. Depending on circumstances, this process can either stimulate or inhibit 
star formation (McKee & Ostriker 1977, Lada 1985, Wolf & Durisen 1987, Spitzer 1990). In a cold, dense cloud, it can 
stimulate star formation by increasing the pressure at the periphery and thus causing collapse in the interior. On the other 
hand, close to the supernova the shock heats the gas and evaporates clouds (McKee & Cowie, 1975). This eventually turns 
off star formation. Supernovae also enrich the gas with metals, which drastically affects cooling rates and the conditions for 
star formation. Much of the energy of supernovae may contribute to the velocity dispersion of cold clouds (McKee & Ostriker, 
1977; Cowie, et al., 1981). 

A proper description of the above phenomena in the context of galaxy formation would not be complete without an 
adequate treatment of gas dynamics. Hydrodynamic simulations in a cosmological context were pioneered by Evrard (1988), 
Cen et al. (1990), and Katz & Gunn (1991), but for a comprehensive discussion of the different numerical methods we refer 
the reader to the comparative paper of Kang et al (1994). The heating of the gas by supernovae and its role in modulation 
of galaxy formation have been studied by a number of authors (Baron & White 1987, Cen & Ostriker, 1992; Katz, 1992; 
Navarro & White 1993; Kauffmann, White & Guiderdoni, 1993; Metzler & Evrard, 1994; Mihos & Hernquist, 1994a; Cole 
et al 1994). Navarro & Steinmetz (1996) studied the effects of photoionization on angular momentum and cooling of disk 
galaxies without including the feedback of energy due to star formation and supernovae and found that photoionization alone 
was not adequate for proper understanding of the " overcooling" phenomenon. Taken together, these previous studies tend to 
confirm the importance of including a realistic description of the star-gas interaction loop including the mechanical effects of 
supernova explosions and multiple gas phases. This point will be addressed further in the conclusions. 
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This paper is organized as follows: In §2 we derive the system of equations describing our multiphase model. In §3, we 
express the equations in a form suitable for application of numerical techniques. For solution of the hydrodynamic equations, 
we have used the piecewise parabolic method (PPM; Colella & Woodward 1984); the collisionless components are treated 
by particle mesh methods (Kates, Kotok, & Klypin, 1990; Klypin & Kates, 1991). The PPM module was written by one of 
us (A.Kh). Several tests of the code are reported in §4, together with parameter studies concerning the effects of supernova 
feedback and enhanced cooling. §5 gives the parameters for a set of CDM simulations and describes algorithms for finding 
galaxies and computing their colors and luminosities using a stellar population synthesis model. Results are discussed in §6, 
including trends in luminous matter vs. total mass, luminosity vs. total mass, colors, and the importance of environmental 
effects. Although the results need to be confirmed by more simulations, we conclude in §7 that the supernova feedback loop 
in combination with hydrodynamics and gravitation lead to a rich variety of dynamical processes in galaxy formation. These 
effects act as "hidden variables" and contribute significantly to the scatter of phenomenological relationships between the 
various observable properties of galaxies. We assume Ho = 50 km s _1 Mpc _1 throughout this paper unless otherwise stated. 



2 PHYSICAL ASSUMPTIONS AND MOTIVATIONS 
2.1 Multi-phase medium 

The matter in the simulated Universe consists of four phases. 1) The dark matter (labeled by a subscript "dm") in th form 
of weakly interacting collisionless particles is the main contribution to the mean density of the universe (Qdm = 1 — fib). 
The baryonic component is described as a medium consisting of three interacting phases: 2) hot gas (labeled by subscript 
h, T h > 2 x 10 4 K), 3) gas in the form of cold dense clouds (subscript c, internal temperature T c = 10 4 K) resulting from 
cooling of the hot gas, and 4) "stars" (subscript *), formed inside cold clouds and treated as collisionless particles. Thus, the 
total density p(r) is the sum of four components: 

P = Pdm + Ph + Pc + P*- (1) 

This picture is consistent with work on models of galaxy formation and evolution (see for example Nulsen 1986; Thomas, 
1988a,b; Hensler and Burkert, 1990; Daines et al 1994; Nulsen & Fabian, 1995) and in this context appears to be superior to 
a treatment of the gas component as a one-phase medium. 



2.2 Cooling of the hot gas 

A proper treatment of cooling requires consideration of processes occurring below the numerical resolution of the code. Despite 
certain simplifications, our picture relies heavily on the comprehensive theory of the interstellar medium of McKee & Ostrikcr 
(1977), who provided a detailed description of the physics in dense regions. This theory presupposes the existence of cold 
clouds forming due to thermal instability in approximate pressure equilibrium with the surrounding hot gas. Heat is conducted 
from the hot to the cold phase through warm interfaces surrounding the clouds, but sufficiently large clouds persist. Despite 
the existence of multiple phases with different temperatures, the rate of cooling can still be computed from knowledge of the 
mean parameters of the hot gas. 

Let us denote the local cooling rate of the plasma due to radiative processes by dE/dt = — A r (p, T), where p and T here 
represent the true local values of the gas density and temperature. According to the theory of McKee and Ostriker, the rate of 
energy loss expressed in terms of average gas density pn and temperature T h will be higher than the nominal rate — A r (p h , T h ) 
by a cooling enhancement factor C: dE/dt = —CA r (p h ,T h ) (McKee & Ostriker used (3 for this parameter). In this paper, we 
take their estimate C = 10. In our numerical code the averaging happens over a resolution cell. This value of parameter C takes 
into account all effects resulting from unresolved density inhomogeneities including radiation from interfaces between the cold 
clouds and the surrounding hot gas. Cooling depends also on the chemical composition of the gas. In order to incorporate the 
effects of metallicity into our code, we assume solar abundance in regions where either thermal instability is present (assuming 
the gas density exceeds an appropriate threshold; see below) or previous star formation has occurred. Otherwise, the region 
is considered not to be enriched by metals, and the cooling rate for a gas of primordial composition (Fall & Rees 1985) is 
assumed. In §2.6, we will describe our procedure for choosing cooling rates in detail. 

The ionization of the intergalactic and interstellar medium by UV photons emitted by quasars, AGNs, and by nonlinear 
structures in the process of formation has two important consequences for our code. First, we assume that the gas outside dense 
regions (p gas < 2(p gas ), where (p ga s) is the mean cosmological gas density) is ionized and has a temperature of T = 10 3,8 K, 
consistent with recent results (Giroux & Shapiro, 1996; Petitjean, Miicket, & Kates, 1996; Miicket et al., 1996). Secondly, in 
low to moderate-density regions, heating due to the ionizing flux suppresses the thermal instability (Miicket & Kates, 1996) 
and thus the formation of cold clouds (and ultimately stars). Hence, we allow cold clouds to form only if p gas /pcr > T> ■ f2b ar , 
with V ~ 50 - 100. 
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In addition to the above processes we include the energy loss due to Compton cooling, given by Acomp = 7x 10 36 njjT e a 4 , 
where uh is the number density of hydrogen ions, T e is the electron temperature, and a is the expansion parameter. 

2.3 Transfer of gas to the "cold" phase 

In this model, there are two mechanisms for gas to leave the hot phase and enter the cold phase. First, if the temperature of 
the hot gas drops below a threshold temperature Tij m = 2 x 10 4 K, and if p ga s/pcr > T> ■ fibar we immediately transfer all the 
hot gas to the cold phase, thus making it available for star formation. This process gives the following terms in the continuity 
equations for the gas: 

(-Tr) = ~T~^-< ^cool = » , > T> • fibar T < Tii m , (2) 

V It / hot^cold tcool Oe^/Ot per 

where t coo \ is the cooling time, and e h is the thermal energy per unit mass. The transfer may be treated as immediate if 
the cooling time is very short compared to a timestep. This typically happens at T ~ 2 x 10 4 K. A second way of transferring 
gas to the cold phase is by sufficiently rapid thermal instability. The temperature range for creation of cold clouds depends in 
reality on the ionizing flux and the local density (Miicket & Kates, 1996). Here, we model these restrictions by checking the 
condition T < Ti ns t, where in this paper Ti nst — 2 x 10 K. To estimate the rate of growth of mass in cold clouds, we assume 
that the energy emitted by the hot gas was actually lost by the gas, which was hot and became cold. (Note that according to 
McKee & Ostriker (1977), the filling factor of the cold clouds is small.) If eh and e c are the internal thermal energies per unit 
mass of hot and cold gas, then the change of energy of the system due to radiative cooling is 

d(p^ + p c e c ) \ = _ CAr(phjThh (3) 

/ cooling 

where C is the enhanced cooling factor due to unresolved dumpiness. We take e c = constant, which means that the cold gas 
cannot cool below Tii m . For an ideal gas, these assumptions imply extra terms in continuity equations: 

/dph\ = _( d fl) CA r (p h ,T h ) 

V dt ) therm. inst \ dt j therm. inst 7 e h — ^c 

where 7 is the ratio of specific heats. 

Because of the frequent exchange of mass between the hot and cold gas phases, it is reasonable to consider the hot gas 
and cold clouds as one fluid with rather complicated chemical reactions going on within it. Thus, we follow the motion of the 
hot component and integrate the change of the total density of the gas (index gas) : 

Pgas = Pc + Ph- (5) 

The cold and hot gas are assumed to be in pressure equilibrium locally. So, P gas = P c = P h = (7 — l)u h , u h — phth, where u 
and e are internal energies per unit volume and per unit mass correspondingly. 

The velocity associated with the fluid involves an averaging by the code over a cell size. We imagine that, in a multiphase 
medium with supernovae, the hot gas is in reality "windy" at scales below the cell size and similarly that the cold clouds 
have a significant velocity dispersion (McKee & Ostriker, 1977; Cowie, et al., 1981; Hensler and Burkert, 1990.) We therefore 
associate an "effective" temperature with the existence of the cloud component which in the present paper is taken simply 
equal to the ambient hot gas temperature, i.e., there is a pressure Pdisp oc p c Th associated with the cold gas component. 
Although this pressure may represent an overestimate, the overestimate in the total pressure (hot plus cold) would only be 
significant if the cold gas fraction could exceed the hot gas fraction. However, in our model this situation never arises due to 
the heating by supernovae as explained below. 

2.4 Star formation 

Star formation, which is assumed to occur only in cold clouds, leads to a decrease in the density of the cold component given 
by 

(if) = "?' (6) 

\ ai / star — formation 

where we take a fixed characteristic star formation time constant t* w 10 8 yr . The actual time-scale for star formation can 
exceed t* and will depend strongly on the rate of the conversion of hot gas to the cold phase. The lifetime of massive stars 
is about 10 7 yr or even shorter; thus most of stars with M* larger than (10 — 20) Mq will explode as supernovae during one 
timestep. Stars in the mass interval from (5-7)M© to 10M Q will explode on a longer time-scale, but they produce less energy, 
and therefore in view of the uncertainties in supernova energy, the energy input due to these stars is not included here. Due 
to the explosion of massive stars, the rate of growth of mass in the form of long-lived stars is decreased by a fraction /3: 
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dp* _ (i-P)Pc (7) 
dt t, ' 1 ' 

The precise value of /3 is sensitive to to the form of the initial mass function (IMF), particularly its low-mass limit. The 
principal source of uncertainty is that this low fraction of mass in massive stars relies on the present IMF, whereas we would 
like to apply it to very early epochs of galaxy formation. The main difference is in the abundance of heavy elements. It is 
possible that in a hydrogen-helium plasma, only supermassive stars could be formed. At the same time, enrichment of the 
medium could proceed very fast, changing the chemical composition of gas and, as the result, the IMF. In view of these 
uncertainties, we assume here that f3 is constant in time, while keeping in mind that other assumptions are permitted and 
might lead to different results. Assuming a Salpeter IMF, we estimate the fraction of mass in stars with mass larger than 
WMq to be /3 = 0.12. 

2.5 Effects of supernovae 

Evaporation of cold clouds is an important effect of supernovae on the interstellar medium (McKee & Ostriker 1977, Lada 
1985). We incorporate this effect by supposing that the total mass of cold gas heated and transferred back to the hot gas 
phase is a factor A higher than that of the supernova itself. We will refer to A below as the supernova feedback parameter. 
Realistic evaporation is a much more complicated process. The rate of evaporation could depend on, among other things, the 
energy of the supernovae, the cloud spectrum, and the ambient density. McKee and Ostriker (1977) give an estimate of the 
evaporated mass which scales with the energy of the supernova E and the gas density as E 6 ^ 5 ■ n^ 4 ^ 5 . However, at present 
we do not include the dependence on the gas density, because here rih refers to the local density, which is relevant for the 
propagation of the supernova shock, whereas our nh is the mean density of the hot gas averaged over quite large volume. 

Assuming that the energy input due to supernovae is proportional to the total mass of supernovae, we obtain for the rate 
of mass exchange due to evaporation 

dph\ _ ( dp c \ _ Af3p c ^ 



dt J evap V dt J cvap t* 

Correspondingly, the net energy supplied to the hot gas phase by supernovae is: 

(^). < 9 > 

Here esN ~ {Esn)/{Msn), where for the supernova energy and mass we take the values (Esn) = 10 erg, (Msn) = 22Mq , 
respectively (Salpeter IMF) . Admissible values of the supernova feedback parameter A are constrained by the condition that 
the energy of a supernova explosion be larger than the energy required for the evaporation of cold clouds. This restriction 
leads to esN > A ■ (eh — £c). For the above supernova energy and mass and for hot gas temperature Th ~ 10 6 K, this condition 
implies the restriction 1 < A < 250. We will discuss implications of varying the supernova feedback parameter A and more 
detailed evidence constraining its value in §^. 

Of course, A could depend both explicitly on time and on local chemical abundances. In the absence of a model for these 
effects, we take constant A, which nevertheless includes the main effects and together with the other effects (energy input due 
to supernovae and radiative cooling) gives a rather complicated nonlinear picture of star formation. (See e.g. Thomas, 1988a; 
Navarro & White 1993 for a different way of modelling the SN feedback). 

2.6 Model equations 

Expressed in expanding coordinates (a is the expansion parameter), the system of equations to be solved numerically thus 
comprises: 

(i) Equations for the dark matter component (collisionless particles): 
dp dx p 

dF = - v ^' Tt = ~^ (10) 

where p = a 2 x and x are the momentum and the position of a dark matter particle, respectively (Peebles 1980). 

(ii) Equations for the "stars" (collisionless particles) with a source term due to the production of new stars in cold clouds: 

dp* dx, p, 

~dT = ~ v< ^ -dt = ^ ' 

^+3^U + ivp^ = iI^, (12) 
at \aj a t* 

(iii) The continuity equation for the gas: 
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(iv) The system of hydrodynamical equations for the "hot" component: 

Ot \aj a U 7£h - e c 

dvh 1 ,__.._/ d \ „ 1 _ „ 1 _ , ,„ 

-^ + -(«h-V)i/h+ (-)vh = VP h --V^, 15 

ot a \aj p h a a 

dE ^ o^U , ivr^/ir , _ -fe h C ■ A r {p h ,T h ) + 



+ 2 (±) Eh + l V{{Eh + Ph){Jh) = _ J^C ■ A r ( Ph ,T h )+ _ 
Ot \aj a je h - e c 



+ ^[e S N+Ae c ], (16) 

where _Eh is the total (thermal plus kinetic) energy of the hot gas per unit comoving volume (E = pe + pv 2 /2), fcTh = 
(7 — l)nnf/i m £h, run is the mass of the hydrogen atom, p, m is the molecular weight per particle, k is Boltzmann's constant, 
and Ph = (7 - l)PhCh 
(v) The Poisson equation: 

A0 = 47rGa 2 (p dm + p» + p gas - p b ), (17) 

where p b is the mean density of the Universe. 

The collisionless dark matter (dm) and "star" (*) components satisfy their continuity equations automatically by virtue 
of the particle-mesh method. Note that the pressure Pdi sp associated with the velocity dispersion of the cold clouds has been 
implicitly included in Eq. (15), since Pdi sp oc p c T h , P h oc phT h , and VPtot/ptot = VP^/p^. As noted earlier, the contribution 
of cold clouds to the density and the pressure turn out to be only a few per cent. 

The above system of equations is solved in two different regimes: If the temperature of a volume element satisfies 
T h > 2 x 10 5 K and if p c is zero, then we suppose that the element has not yet participated in star formation. In this case the 
constant C is set to unity, and primordial chemical composition is assumed. If either the temperature is lower or p c is not 
zero, or there is already a significant star density present, then solar abundance is assumed with the cooling rate given by 
Raymond, Cox, and Smith (1976). Compton cooling is always included. We plan to make the cooling rates more realistic by 
adding to the model the density of heavy elements produced by supernovae and by modeling the effects of photoionization in 
a self-consistent way as discussed in the conclusions. 



3 NUMERICAL TECHNIQUES 

3.1 Equations in expanding rescaled coordinates 

In order to apply the PPM method (Colella & Woodward 1984) for gas dynamics in a cosmological framework, we seek 
transformations from the above physical quantities to rescaled quantities satisfying equations that have the "usual" form 
of equations of gas dynamics (i.e., in the absence of an expanding background). We adapted the transformation used by 
Shandarin (1980). First, we change to expanding coordinates. Auxiliary dependent variables (denoted by a tilde) are related 
to the physical variables as follows. Density: p — a 3 p; peculiar velocity: v = av; pressure P — a 5 P; internal energy per unit 
volume u = a 5 u; internal energy per unit mass e = a 2 e. Second, we introduce a new time variable 

b{t) = ^-(l-a- l ' 2 {t)), (18) 

where Ho is the Hubble constant now. For convenience, b(t) was chosen to vanish at redshift zero (a = 1). Thus, b is negative 
for a < 1. (In a nonfiat background, the formula for b would be changed; the definition of b would be da = aa 2 db.) 
In terms of the new rescaled variables, the Poisson equation and the Euler equation take the form: 

2 - A-kG 

V (f> = — — (pgas + /5dm + /5, - Pb), (19) 

^ + (4-V)4 = -|vA-a 2 Vf (20) 
Ob p h 

The left-hand side of the continuity equation is transformed back to the usual form (i.e., no terms containing a), but the 
terms on the right-hand side are multiplied by a factor a 5 . The same thing happens to the energy equation: the left-hand side 
takes the usual form, while the terms on the right-hand side are multiplied by a power of the expansion parameter (a 7 ). In 
rescaled variables, the energy equation and continuity equations for the hot gas, for the gas component, and for "stars" take 
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the form: 



dE b 
db 



+ V ((E h + P h K) 




(21) 



-?Jb +VphVh 



2PA Pc aCpiL(T h )/{fi H m H ) 2 



(22) 



U 7?h - £c 



dp. 



db 



+ VpgasCh 



3 2 (i - p)h 



(23) 




2 (1 ~ P)Pc 

a 

u 



(24) 
(25) 



Here we have defined A r (ph,Th) = (ph/ Hhtuh) 2 L(Th), where p^r is the molecular weight per hydrogen atom and inn is the 
mass of the hydrogen atom, and Dcomp = 7 x 1CT 36 / (fiHrnn)- 

3.2 Solution of the hydrodynamical equations 

To calculate the flow of the ordinary (collisional) matter, described here as a multi-component medium, we apply a hydrody- 
namical code based on the PPM technique. One of the main advantages of this Eulerian code is that it does not involve an 
artificial viscosity, but still treats shocks very accurately (for example, it does not generate post-shock fluctuations, and shocks 
are only one zone thick). The method is very fast and highly parallelizable. This algorithm is applied to our gas-dynamic 
equations in one dimension at a time. Multiple dimensions are treated by directional timestep splitting, whereas local processes 
(heating, cooling, star formation, etc.), which involve various components of collisional matter, as well as self-gravitation and 
gravitational interaction with dark matter, are treated using process timestep splitting (Oran & Boris 1986). Composition 
variables (such as the concentration of heavy elements or the density of the gas p g as) are advected through the mesh by the 
same PPM algorithm. 

The time-scale of the cooling and star-formation processes usually is shorter than the dynamical time-scale. In our code 
the time step for advection of the gas and motion of the dark matter is equal for all zones, but cooling and star-formation 
processes are integrated with a timestep adapted for each zone. 

4 TESTS AND PARAMETER STUDIES 

4.1 Hydrodynamic tests 

In order to check the accuracy and efficiency of the present implementation, we first verified proper behavior in a number of 
classical tests, two of which are presented here: Fig. 1 shows the density and pressure for the Sod (1978) problem for the decay 
of a discontinuity with the parameters p L = 1, P L = 1, = 0, Pr = 0.125, Pr, = 0.1, Mr = 0, 7 = 1.4. All the important 
features of the flow (the rarefaction wave on the left, the contact discontinuity in the middle, the shock on the right) are well 
represented. Fluctuations of physical parameters in the vicinity of the shock and contact discontinuities do not occur. The 
shock and the contact discontinuity are spread over only one cell. These are important desirable properties of the code. 

In Fig. 2, we present the results of a strong "2-D" cylindrical explosion in an ideal gas with 7 = 5/3 as performed on 
100x100 mesh with periodic boundary conditions. The explosion energy is initially deposited inside a small cylinder of radius 
2.5 cell widths located coaxially to the 2-axis with P= 10 4 . The pressure outside the small cylinder is initially 0.1. The initial 
density is p — 1 everywhere. The lower panels show the density at different moments. In particular, they demonstrate that 
oblique shocks are stable. In the left top panel, the density of every cell in the simulation is plotted against its distance from 
the center of the explosion. The narrowness of the curves demonstrates that the shock propagates in all directions at the same 
rate, and that the thickness of the shock stays « 1 cell even after significant evolution. The top right panel shows how the 
shock expands. Note that in the case of cylindrical symmetry the appropriate dimensionless variable is £ = r x (a/(Et 2 ))~ 1 '' 4 , 
where a is the surface density, rather than the usual "Sedov" -dimensionless variable £ = r x (p/ Et 2 )^ 1 '' 5 . The resulting 
expansion rate is, according to the Sedov solution, r oc t 1 / 2 . 

4.2 Supernova feedback and cooling enhancement parameters 

The evolution of the different baryonic components provides information useful in understanding the role of the supernova 
feedback parameter A and the cooling enhancement factor C as control parameters in our nonlinear star-gas model. Cooling 
and star formation were studied in a single cell, i.e., with gravitation and mass flow switched off. Fig. 3 shows the time 
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evolution of the relative mass fraction of the different phases for models with varying supernova feedback parameters A = 1 
(no evaporation), A = 20 (w 10 per cent of supernovae energy goes to evaporation) and A = 100 (« 50 per cent goes 
to evaporation). The gas was assumed to have solar composition and initial temperature T = 10 6 K and density hh = 
6 x 10~ 4 cm~ 3 , typical of a galaxy halo. The system was assumed to be in the regime of thermal instability with cooling 
enhancement factor C = 10, which as explained above accounts for unresolved substructure in the interstellar medium. 
Independent of A, the hot gas initially cools with a characteristic cooling-time of (2 — 3) x 10 7 yr . At t ~ 5 x 10 6 yr , the first 
stars form and begin to produce supernovae. The supernovae evaporate cold clouds, which is evident for example in the case 
A = 20 as a rapid drop in the density of the cold gas and as a jump in the density of the hot gas. Although at early epochs 
(t < 10 8 yr) larger feedback parameter A results in a smaller fraction of mass in cold clouds and in "stars", the trend later 
reverses: high A ultimately causes more of the mass to end up in stars during the time available, because the temperature of 
the hot gas stays lower and it cools faster. The star-formation e-folding time was 4 x 10 8 yr for A — 100, 10 9 yr for A = 20 
and 3 x 10 9 yr for A = 1. Thus, due to the nonlinear feedback the true star formation time scale is greater than t„ = 10 8 yr . 

The results of this and other parameter studies illustrate that the supernova feedback parameter A is directly related to 
the efficiency of star formation. Indeed, in the context of a full simulation, the effect of A may be even stronger than implied 
by these studies: Relatively larger values of A imply that a larger fraction of the energy supplied by supernovae will convert 
gas from the cold to the hot phase. As seen in Fig. 3, on the short term this conversion will temporarily inhibit star formation, 
but by causing a smaller pressure gradient it will allow at least the larger galaxies to keep most of their gas. On the other 
hand, relatively small values of A imply that a large fraction of energy supplied to the gas by supernovae will heat the hot 
phase to even higher temperatures locally. The ensuing pressure gradient will enhance the expulsion of gas, even from large 
galaxies, and retard star formation for a time comparable to a dynamical time scale. It is also clear that the dynamics of 
star formation will be sensitive to the effects of environment such as mergers, gas flow, tidal forces, etc., because these effects 
strongly influence the amount and temperature of gas available for star formation. 

The value of A should be chosen to allow for such phenomena as "young, active galaxies" and the transport of metal-rich, 
hot gas from galaxies to the intergalactic medium. Because we are assuming A constant in time, our choice will represent a 
compromise between the necessity of allowing for really active regions of star formation (low values of A) and low efficiency 
of conversion of gas to stars inside cold clouds (high A). These constraints suggest a value for A between 20 and 200. Another 
constraint comes from the condition that the total mass in luminous matter ( "stars" ) be comparable with the observed value 
of l%-2% of the critical density of the Universe. From our simulations we infer that choosing the feedback parameter in the 
range A — 100 — 200 gives about the correct amount of mass in "stars" . 



5 SIMULATION OF THE EFFECTS OF FEEDBACK MECHANISMS ON GALAXY FORMATION 
5.1 Cosmological model and initial conditions 

The initial fluctuations were computed as a realization of an expanding flat CDM model (according to the approximation of 
Bardeen et al. 1986) normalized to the rms mass fluctuation in a sphere of radius 8/i _1 Mpc to be as = 0.67 as estimated by 
the linear theory. The ratio of baryonic to total density was taken to be f^j, = 0.08. 

We present results of ten numerical simulations whose parameters are given in Table 1. We used 128 3 particles and 128 3 — 
256 3 cells in the simulations reported here. The simulations were carried out on SGI Power Challenge parallel supercomputers 
at NCSA and at the Centro Europeo de Parallelismo de Barcelona (CEPBA). About 700-1100 timesteps are needed to evolve 
a simulation to z — with timesteps chosen to be 0.6 of the maximum according to the Courant-Friedrich-Levy condition. 
Due to the high degree of parallelization of the code, one time step requires only about one minute of real time on the SGI 
Power Challenge with 4 90MHz R8000 processors; it takes about 100 node-hours to run a simulation of 1000 steps. For a 256 3 
simulation, one time step takes about 8 times as long, and 1 GByte of main memory is required. 

The simulations were performed for a variety of different parameters and resolutions. The three last lines of Table 1 give 
the fraction of mass Oi um i nous converted to stars and the redshift Zi/ 2 at which f2i um i n ous was half of its value at redshift zero. 
(We assume Ho = 50km s" 1 Mpc _1 throughout this paper unless otherwise stated.) 

The choice of box size clearly involves a compromise among different considerations. Sufficiently small cell size was the 
overriding consideration for this paper. The box size of 5 Mpc actually chosen for Models 1-5 implies a cell size of 39 kpc, 
which is just large enough to resolve galactic scales of interest. At the same time, a box of this size contains quite a few 
galaxies of different masses at z = 0. It also typically contains a filament and thus incorporates at least this aspect of the 
observed large-scale structure. Models 1,3,4, and 5 were run with exactly the same initial conditions. Model 2 had a different 
realization of initial fluctuations, but otherwise the same set of parameters as Model 1. Model 10 was the same realization as 
Model 2 (same box size and phases), but with 8 times as many cells (size 19.5 kpc). 

The smaller the box size, the greater the influence of the boundary conditions: The size of the largest structures in the 
box should be small compared to the box, because mirror images due to the cyclically symmetric Green's function would 
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otherwise distort the results. Also, a smaller volume clearly reduces the number of galaxies to be expected at the end of the 
simulation. Nevertheless, in order to check the sensitivity of our results to the effects of resolution, we ran another series of 
models (Models 6-9) with comoving cell size 11.7 kpc (box 1.5 Mpc). 

Models 1-6 were evolved up to z = 0, while Models 7-10 were evolved until z — 1. 



5.2 Simulation parameters and constants 

We now discuss and summarize the specific parameter values used here. 

The supernova mass fraction j3 defined in §2.5 was computed assuming the Salpeter initial mass function dN* /d(log M) oc 
M -1,35 for masses between O.IMq and IOOM0 . Assuming that stars with M > WMq explode as supernovae, we obtain 
P = 0.12, the value adopted here. One can consider the effect of steepening the IMF. For example, an exponent of -1.6 would 
lead to /? = 0.05. Note that at the low- mass end, there is evidence for a decrease in the slope of the IMF (Miller & Scalo 1979). 
Thus, taking the lower mass limit to 0.5Mq , one would obtain (5 = 0.13. 

In order to arrive at a value of the characteristic star formation time constant introduced above i» = 10 s yr, we note that 
i) even massive stars live for tyis = 3 x lO r [M,/lOAf0 ] _1B yr ~ (1 — 3) x 10 7 yr before they explode as supernovae, and ii) 
star-formation occurs in different clouds; it is reasonable to suppose that the production of stars inside different clouds in one 
simulation cell is not synchronous. 

In introducing a single constant t», it is important to recall that in our picture, star formation is regulated by the 
availability of cold gas clouds. In codes that do not provide an independent estimate of the cold gas component, a reasonable 
procedure is to estimate a local dynamical or cooling time of the gas component in order to infer the amount of cold gas 
available (e.g. Katz, 1992; Navarro & White 1993). In our code, the dependence of the star formation rate on the cooling time 
is already implicit in the dynamics of the multiphase model. 

In §^|, we saw that the star formation rate depends strongly on the supernova feedback parameter A. In Models 1 -7 and 
10, we take A = 200, corresponding to very efficient star formation and for comparison A — 100 in Models 8-9. 

"Stars", represented by particles, were assigned a mass implied by Eq. (jj). Once produced, a "star" particle moves 
like a collisionless particle. (In fact, these particles are best thought of as representing starbursts.) As explained in §2.2, 
photoionization will tend to prevent cooling in low-density regions. To incorporate this effect, a threshold of T> = 50 — 100 
times the background baryon density (see Table 1) was imposed for converting hot gas to cold clouds in a cell (i.e., switching 
on star formation in a cell). This prescription also has the advantage of reducing the computational cost of excessive low- mass 
starbursts. 



As explained in §2.6, solar abundances were assumed if either star-formation conditions were met or stars were already 
present in a cell. Now, the presence of formed stars in a volume element is an important indicator that the gas was enriched 
by metals. However, in simulations reported earlier (Klypin, Kates, & Khokhlov, 1992; Yepes et al., 1995 , 1996a), the cooling 
regime was determined from star-formation criteria alone. Note that in these earlier simulations the hydrodynamic module 
was based on the method of flux corrected transport (Boris & Book, 1973, 1976). Hence, as an additional control, it was 
useful to run a few models (5,8, and 9, marked by asterisks) without taking into account the presence of star particles when 
choosing the cooling regime (primordial/solar composition). These runs also illustrate the sensitivity of results to the chemical 
composition. It turns out that the effect of taking primordial composition instead of solar hardly affects the early evolution 
of Models 5, 8, and 9 compared to the others: Since the temperature of the gas is still comparatively low, the gas cools 
fast with either composition, producing cold clouds. It also has little effect in small galaxies at any epoch, because in small 
galaxies, winds (pressure gradients) expel the heated gas in either case. But the evolution at later times (z < 1) is different. 
The temperature of the gas in large galaxies grows with time because of the energy released by supernovae and because more 
massive galaxies are capable of holding on to the hot gas. Thus, the cooling time gradually increases in massive galaxies. Now 
the system becomes sensitive to the assumed chemical composition. Eventually, the density of cold gas will drop below the 
threshold - a few per cent of the total gas density - needed for a star formation to occur. At this moment, the models in which 
the local presence of stars is not taken into account switch to primordial cooling, which effectively shuts off star formation for 
the galaxy after z = 0.5. If the local presence of stars is taken into account, the code is more likely to switch to solar cooling 
and hence higher cooling rates. Thus, the galaxy tends to produce cold clouds and new stars at a steady rate. In spite of this 
qualitative difference, which dramatically affects star formation rates and thus the colors of large galaxies, the overall effect 
on the total mass in stars and on the epoch of star formation Zx/2 was still not that significant (compare models 4 and 5 in 
Table 1). 



5.3 Galaxy identification algorithm 

Since the code produces separate estimates of all four density components (dark matter, hot gas, cold gas, and stars), the 
information contained in these four density fields could in principle be used in a variety of ways to define "galaxies" . From the 
point of view of comparison with observation, one is most interested in the distribution of visible matter, suggesting that the 
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star density should be viewed as the most important quantity to be analyzed in finding galaxies. Nevertheless, identification 
of dark matter halos proved to be the most useful procedure for the purposes of this study: In order to study the relationship 
between galactic morphology (e.g., luminosity, spectral properties) and location or environment, we obviously need to include 
low luminosity but relatively massive halos in the analysis. 

For this reason, we begin by constructing the total density p on the grid and searching for local maxima in p such that 
Sp = p — pb exceeds a fixed threshold density pthresh- (Here we took pthrcsh = lOOpcrit.) Next, we place spheres of r sp = 2 cell 
radii around maxima and determine the center of mass (CM) of the particle distribution (dark matter and stars) within the 
spheres. (In view of limited resolution and the low mass fraction of gas, inclusion of the gas density would have little effect.) 
We then translate the center of the sphere to the CM and compute the mass M sp and the overdensity (8p/p) sp within the 
sphere. If (Sp/p) sp > 200, we add the halo to our catalog. If (8p/p) sp < 200, we reduce the radius to r sp = 1 cell radius. If 
(8p/p) sp > 200 in this sphere, then we add the halo to our catalog; otherwise we do not include the halo. 

For models with cell size 39 kpc, the lower galaxy mass limit for r sp = 2 is roug hly 2 x 10 10 M Q , while for r sp = 1 the 
limit is roughly 2.5 x 10 9 M Q . We have checked that no overlapping halos are identified with this algorithm. 

5.4 Stellar Population Synthesis 

We wish to estimate the spectral energy distribution (SED) S(\, t) at an arbitrary epoch for each halo (referred to here as 
"galaxy") in the catalog. The SED will then allow us to compute luminosities and colors, to be analyzed below. For each 
galaxy, we have a list of particles representing the stellar component. Associated with each particle is the mass and time of 
formation of the "starburst." Consistent with our assumptions of constant IMF and star formation under conditions of solar 
chemical abundances, these two parameters allow us to compute the evolution of the SED from a stellar population synthesis 
model. 

A number of stellar population synthesis models have been proposed (e.g. Tinsley 1972, Bruzual 1983, Guiderdoni & 
Rocca-Volmerange 1987, 1990, Buzzoni 1989, Chariot & Bruzual 1991, Bruzual & Chariot 1993). Among them, we have 
chosen the model of Bruzual & Chariot (1993) that describes the time evolution (between and 20 Gyr) of SED's for a 
burst of star formation under conditions of solar metallicity and a Salpeter IMF, in accordance with our model. It covers 
the relevant stages of stellar evolution, including the asymptotic giant branch (AGB) and post-AGB stars, and it is based on 
up-to-date stellar evolution calculations. 

For each galaxy in the catalog, we may compute S(\, t) from 



where 3>(Tj) is the mass of stars in the halo produced at timestep n and T{\,t) is the SED of a starburst of 1 Mq after 
a time t. Then we simply convolve S(\,t) with the filter response function Rf(X) to obtain the absolute luminosity Lf(t) 
in the given band. Combining the Lf(t), we then obtain the evolution of the color index of the galaxy. Of course, the color 
of a galaxy that would actually be measured will be influenced by other factors such as the interaction of starlight with the 
surrounding plasma. 



6 RESULTS 

We begin our discussion of the results by comparing some global statistics concerning observational quantities with estimates 
obtained from galaxy catalogs. For example, in Table 2 we give the average B-band luminosity density for Models 1 — 5(5 Mpc 
box). These values were obtained by integrating all the light from the stars at z = 0. Depending on the different parameters, 
the luminosity density ranges from 7 — 11 x 10 7 Lq/Mpc 3 . Observational estimates from the APM redshift sample (Loveday et 
al., 1992), which covers the largest volume surveyed to date, give a B-band luminosity density of 6.2 x 1O 7 L /Mpc 3 . Taking 
into account that there could be systematic photometric errors in their estimates (e.g. Ellis et al. 1996; Bertin & Dennefeld, 
1996), a factor of 2 increase cannot be ruled out. This would imply that the amount of light produced in our models is not 
far from what is estimated for the Universe. Bolometric mass-to-light ratios M/L can also be computed from our simulations. 
They range from 24 to 85 Mq/Lq when averaged over a volume of 5 Mpc on a side. On the other hand, the star formation 
rate at low redshifts in our simulations was a factor of « 4.5 higher than the most recent estimate of O.O13M /yr/Mpc 3 
(Gallego et al. 1995) for our local Universe. 

An important characteristic can be deduced from Table 2. Although as we shall see below the observational properties 
of individual halos can be strongly affected by different values of the "chemical" parameters, the global statistics of the 
system are less sensitive to these values. Figs. 4 and 5 compare Models 1 and 2, respectively, which have the same chemical 
parameters but different realizations of the initial fluctuation spectrum. In Figs. 4a and 5a, 3D views of the dark matter 
distribution in the respective simulation boxes at redshift z — are shown. In Figs. 4b and 5b, the distribution of "star" 
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particles is illustrated. (Note however that the figures do not take into account the different masses assigned to different 
"star" particles.) Comparison of the figures reveals that there are many dark halos, especially smaller halos, which have not 
succeeded in producing stars. "Biasing" is suggested by visual inspection of these figures, in that the more massive objects 
contain a disproportionate number of star particles. The star concentrations trace all of the large dark matter halos, but there 
is more variation in the luminosity of smaller-mass halos. There is a striking visual impression that virtually all halos located 
in voids are completely dark. 

Figs. 6a (Model 1), 6b (Model 6), and 6c (Model 10) show contour plots of the temperature distribution T gas as well 
as the densities of luminous matter ("stars"), dark matter pdm, and gas p gas in a thin slice of thickness 1 cell width passing 
through the centers of the most massive galaxies in the simulations. The filamentary structure is again prominent in the gas 
distribution. The most massive galaxy in Fig. 6a (just above the center) is rather large - it has mass 7 x 10 Mq and 3D 
rms velocity 250 km/s. The internal gas temperature is 2 x 10 6 K. The galaxy is surrounded by a large (1/2 Mpc radius), 
low-density "cavern" of hot ~ 10 6 K gas. The cavern is produced by stellar winds blowing from the galaxy and by a converging 
flow of colder gas falling onto that galaxy, primarily from the filaments. The hot gas tends to expand into low-density regions. 

There are two small satellites of the galaxy inside the cavern. In spite of their small mass, these satellites have managed 
to produce a significant amount of luminous matter as clearly seen in the upper right panel of Fig 6a. Compare those small 
galaxies with their counterparts further from the large galaxy (e.g. at x = 3.5, y — 2, or at x = 2, y = 4.7), which do not 
produce "stars". This is an example of biasing depending on the environment. 

Galaxies in Fig. 6b (Model 6 at z — 1) are significantly smaller than the large galaxy in Fig. 6a. The largest galaxy in Fig. 
6b (at the center) has 3D rms velocity 80 km/s. As the result, the size of the hot gas cavern and the temperature of its gas 
are significantly lower. Nevertheless, the basic pattern is similar: star formation predominates in large galaxies and in smaller 
galaxies which are inside relatively dense filaments. Small dark halos outside the filaments (and far from a large galaxy) do 
not generate "stars" and luminosity. 

In Fig. 6c, which represents the high resolution Model 10 at z ~ 1, the pattern of preferential formation of stars along 
the filament continues, even for dark matter halos with a significant amount of gas. The change in resolution does not lead to 
a qualitative change in the picture. 

There is a strong statistical correlation between the density of the luminous matter and the density of the dark matter. 
In Fig. 7a the density of the luminous matter and the density of the dark matter are plotted for every cell in the simulations. 
Points with density contrast less than 200 represent either the periphery of galaxies or lie entirely outside galaxies. They show 
very little correlation of the luminous matter with the dark matter. But high density and high luminosity density regions 
demonstrate quite a remarkable correlation, which we estimate as 

Plum = (pdm/133) 1 ' 7 , (27) 

for z = and approximately 3 times smaller amplitude for 2 = 1. The only models which had significant deviations from this 
correlation are Models 8 and 9 (the latter is not shown in Fig. 7a). Those models were run for A = 100 - one-half of the 
value in the other models. A smaller value of the feedback parameter resulted in significant suppression of star formation in 
high-density regions. Those models had ~ 3 times smaller total amount of luminous mass (Table 1). It is interesting to note 
that changing the cooling rates by ten times (Model 1 vs. Model 3; Model 8 vs. Model 9) does not have much of an effect. 
The scatter of this relation is still significant - a factor of 2-3. It leaves room for the non-local biases observed in Figs. 5 and 
6. 

In Fig. 7b (Model 10), we observe (up to a constant of proportionality) the same trend at redshifts 4, 2.3, and 0.9 as in 
Eq. (I27]). There is no effect of resolution. This provides evidence that the relation could retain its validity at higher densities. 



In Figs. 8,9, and 10 we plot different characteristics of "galaxies" identified with our algorithm (see §5.3): the ratio of 
luminous to total mass M s tars/Aftot as a function of the total mass of the galaxies (Fig. 8), the absolute visual magnitude of 
the halos versus total mass (Fig. 9), and the UBV color-color diagram (Fig. 10). The different plotting symbols in the figures 
indicate the distance to the most massive galaxy. The dashed line in Fig. 8 shows the average baryon fraction fi_g = 0.08. 
For high- mass "galaxies" in Fig.8 (M tot > 3 x 10 10 M© ), there is a trend in the mass fraction in luminous matter as a 
function of the total mass, which we approximate as M sta rs/Aftot ~ 0.035-y/Aftot/lO 11 Mq . The trend indicates that larger 
galaxies are more efficient in converting gas into stars. It may be understood by imagining that the efficiency of gas ejection 
due to supernovae is a strong function of the mass of the galaxy, at least in some range of galaxy masses, because positive 
pressure gradients, which are not primarily dependent on the mass of the object, will be balanced by relatively large inward 
gravitational forces, which are related to the mass of the object. At the very high-mass end, for A = 200, gravitational forces 
already suffice to avoid gas loss, and therefore further increases in mass would have a less dramatic effect on star formation, 
leading to a flattening of the curve. For lower values of the feedback parameter A (Models 8 and 9), the pressure gradients 
are higher because more energy from supernovae goes into reheating the hot gas. The net result is a global decrease in the 
amount of stars formed for some given total mass. This effect is clearly seen in Fig. 7a for Model 8 (lower right panel). (See 
also Yepes et al. 1995, 1996) 

As the mass decreases, the M s tars/A4tot curve first flattens off and then develops a large scatter at low masses. The 
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dramatic variation in the fraction of luminous matter down to zero is clearly an indication that the gravitational potential 
of these galaxies was not deep enough to bind the gas heated by the first supernovae which occurred at early times when 
the galaxies first collapsed. These are the galaxies that did not sustain even the first burst of star formation. Those low-mass 
galaxies with relatively high star mass fractions and/or higher luminosities had their gas replenished by the environment. 

Slightly larger galaxies (masses around 10 10 M© ) had another history. We have traced the star formation rates of few of 
these galaxies (Figs. 11 and 12). These galaxies were able to keep their hot gas for some time: the Universe was denser at high 
redshifts, and the gas was still cooling rather fast and falling back onto the galaxies. But as time went on, the gas was getting 
farther away and its cooling time was getting longer. This resulted in a significant decline in star formation rates at small 
redshifts. Because the star formation rate peaked at high redshifts (z = 3 — 5) for small galaxies, those galaxies are relatively 
less luminous and redder at present. Because of that the flattening in the ratio M BtSlTB /M tot at M ~ 10 10 A/q (Fig. 8) is barely 
seen in the luminosity My (Fig. 9). It is possible to check the effect of spatial resolution on the peak of star formation by 
considering the higher resolution Models 6-10. We find the same trend: the peak of the star formation rate tends to higher 
redshifts for progressively less massive galaxies. 

The dependence of the absolute visual magnitude My on the total mass can be approximated as 

Mv = -18.5-4lQgf-^-V (28) 



1O 11 M 

This relation fits all our models for Aftot > 2 x 1O 1O M0 and My < —15 and has very small deviations of less than 1/2 
magnitude. 

A large spread in the fraction of luminous matter is apparent for the low-mass halos. Figs. 6 give evidence that low-mass 
halos located along filaments are more likely to have a higher mass fraction of stars than those not located on or near the 
filament. Such a correlation would be expected in a scenario in which gas traces the local filamentary structure and is thus 
preferentially available for star formation in those galaxies located along the filament. As explained earlier, for the value 
of the supernova feedback parameter A = 200 used in most of the models, the expulsion of gas from dense regions due to 
supernovae is relatively inefficient compared to evaporation of cold clouds, increasing the amount of hot gas at the expense 
of its temperature. In massive galaxies, this gas remains in the potential well, cools, and falls back into the central regions 
to allow later epochs of star formation. Now, it is true that even for A — 200 supernovae would eject most of the gas out 
of the potential well of a small galaxy, but there is a relatively large reservoir along the filament including cold gas that can 
eventually lead to some future star formation. 

Fig. 10 is a UBV color diagram for galaxies brighter than My < — 15 in our catalog (in the rest frame of the galaxies). 
In each panel, the horizontal and vertical axes represent the Johnson B — V and U — B color indices, respectively, computed 



as explained in §5.4. The long dashed curve represents the UBV color diagram for main sequence stars of luminosity class V. 
For reference, the short dashed line represents fully corrected colors of observed galaxies in the range from Hubble types +10 
(irregular) at the upper left end to —6 (E0) at the bottom right end (de Vaucouleurs, 1977). The typical deviation is about 
0.1 magnitude. 

The simulated galaxies are roughly consistent with the observed relation, although there is a slight blue trend. The blue 
trend is reduced in Model 4 (lower right) compared to Model 3 (lower left). In particular, note that in the lower left panel 
B — V for the most massive galaxy in the simulation is shifted by about 0.2 toward bluer colors compared to the lower 
right panel. The bluer color is associated with a higher rate of star formation at recent times as discussed for Fig. 12 below. 
The color differences may be attributable to the lower density threshold T> of Model 4 for allowing star formation. We recall 
that the density threshold is intended to incorporate in a simple way the effects of heating in lower-density regions due to 
photoionization. 

Note that, according to the Bruzual & Chariot (1993) model used here, a B — V color index greater than 1 as in elliptical 
galaxies would require more than 13 Gyr, and in view of Ho = 50 km s _1 Mpc _1 this possibility is therefore excluded. 

In Fig. 11, the star formation rate (SFR) is plotted as a function of redshift for three galaxies in Models 1 and 6: Star 
formation begins very early (z > 9) in the largest galaxy (solid line), and the SFR rises up to a maximum of ~ 15M0 /yr at 
about 12 Gyrs ago (z = 2), after which it falls steadily towards a present SFR of (2 — 3)M© /yr. In the third-largest galaxy 
(dashed line), star formation begins later, rising to a maximum at about z = 2.5 and then falling off more sharply. Note that 
Model 6 has more than 3 times better resolution than model 1 (5.8 kpc at z = 1) and yet it basically shows the same results: 
the peak in star formation is higher and it occurs later for more massive galaxies. The decline is stronger for smaller galaxies, 
which means that dwarf galaxies (like those least massive shown in this plot) are older and redder than more massive galaxies. 
The periodic bursts of star formation exhibited by the small galaxy (dotted line) are rather typical for the dwarf galaxies in 
our simulations. This behavior is due to the shallow potential wells of these halos, which are unable to retain the gas heated 
by supernova explosions (Yepes et al 1995). The positive pressure gradient drives the hot gas out of the central region, but 
some of it will still manage to remain in the outer parts of the potential well. The gas may then cool, most efficiently in the 
central regions, and subsequently recollapse, giving rise to a new burst of star formation. In this sense, our simulations seem 
to lend support to some previous analytical studies of the nature of dwarf galaxies in the Universe, in which the periodic 
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bursts of star formation for these objects and their dependence on environment were predicted (Dekel & Silk 1986; Babul & 
Rees 1992; Campos- Aguilar, Moles & Masegosa 1993) 

Fig. 12 illustrates the dependence of the SFR on model parameters while holding the phase relationships in the initial 
realization of the perturbation spectrum constant. In the upper panel, the upper curves are for the most massive galaxy, while 
the lower curves are for the least massive luminous galaxy. The solid lines correspond to Model 1 (also seen in Fig. 11); the 
dotted lines are for Model 3 (no enhanced cooling); while the dashed line is for Model 4 (lower star-formation threshold). The 
lower panel is for the high-resolution simulations: the solid lines are for Model 6, the dotted lines for Model 7 (no enhanced 
cooling), and the dashed lines are for Model 8 (lower supernova feedback parameter). At high redshifts (z > 4) galaxies of 
a given mass category behave almost in the same way independent of our simulation parameters. At redshift one, Model 8 
(smaller feedback parameter A = 100) gives a factor of at least five lower star formation rates than the models with A — 200. 
As pointed out earlier, for A = 100 a smaller fraction of the supernova energy is used to transfer gas from the cold to the hot 
component, and the net effect is to raise the temperature of the hot gas more, thus increasing outward pressure gradients. 
Strangely enough, Models 3 and 7 (no enhanced cooling) actually have a somewhat higher rate of star formation at recent 
times as discussed in Fig. 10. All in all, the effect of changing C by a factor of 10 is not as dramatic as one might have 
suspected. 

Above, we noted that the global star formation rate at recent times seems to be higher than the observed rate even 
though the luminosity comes closer to agreeing. A separate but related issue is the even higher star formation rate early in 
the simulations, which reflects the prevalence of early objects due to high small-scale power in the standard cold dark matter 
spectrum. This is also reflected in the very red colors of most of our old dwarf galaxies. 



7 CONCLUSIONS 

Confrontation of scenarios for large-scale structure formation with observations require an assignment of observational prop- 
erties such as luminosity and colors to mass concentrations. The procedure utilized here provides such an assignment using a 
model based on astrophysical processes in the interstellar plasma which (with reasonable confidence) we believe must occur, 
even if the details are not completely understood. While improvements in numerical resolution are always welcomed, at some 
point in any approach it will be necessary to estimate the effects of processes occurring below the resolution on the basis 
of what is known about quantities averaged over the resolution scale. These estimates demand careful consideration of the 
small-scale nongravitational physics, which plays a major role in the most important observable consequences. 

All of these effects contribute to and complicate the relation between the luminosity and colors of galaxies on total mass. 
At sufficiently small scales, the aforementioned strong coupling of morphology to dynamics and thus to physical location poses 
severe problems for any simple-minded biasing prescription such as those frequently used in large-scale structure simulations. 
This conclusion generally supports our working hypothesis that improvements in resolution resulting from improved computing 
resources will only lead to better predictions if the most important processes occurring at small scales are properly modelled. 
As we saw above, the colors of the brightest galaxies found here are dominated by young stars, tending to younger stellar 
types with increasing luminosity. However, it is too early for any sweeping statement. 

The results of this paper underline the rather complicated interaction of gas dynamics, supernova feedback and gravitation 
in the evolution of observable properties of galaxies. Gravitational condensations are of course a prerequisite in order that an 
adequate supply of cold gas be available for star formation. Subsequently, gravitation itself plays more of a passive role: the 
most critical processes regulating the dynamics of star formation (and ultimately the observable properties of galaxies) are the 
various star-gas interactions and feedback loops. The most important of these is the supernova feedback loop, which regulates 
star formation in three different ways, two inhibitory and one stimulatory: First the evaporation of cold clouds reduces the 
efficiency of star formation. Second, the enormous transfer of heat to the hot gas component results in pressure gradients, 
which work against gravitational forces and expel gas, especially from shallow potential wells. This effect strongly couples 
subsequent star formation to hydrodynamics, and in particular provides a possible mechanism for a strong dependence of 
morphology on position. Third, metal production due to supernovae further enhances cooling. This third effect is only seen 
locally in our code, because composition gradients are not advected, but it should in principle cause some additional further 
coupling of the hydrodynamic and supernova feedback loops. 

A relation between the absolute visual magnitude My and the total mass was given in © for Af to t > 2 x 1O 1O M . 
Some of the "scatter" in this relation can presumably be attributed to the environment, but much more work will be needed 
to establish a pattern. 

In this paper, effects of photoionization were included only in a phenomenological way. We would expect photoionization 
to inhibit star formation most effectively at the low-mass end of the mass luminosity relation. There would be several important 
effects. First, one might expect relatively more young stars and therefore increased visual magnitude now. Second, the colors 
would be shifted toward the blue. For larger galaxies, we would expect fewer modifications of the present results. However, it 
seems unlikely that proper inclusion of these effects would cancel out the most important qualitative conclusions of this paper, 
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namely that there is a strong coupling between gas dynamics and the supernova feedback loop, which in turn act as a kind 
of "hidden variables" and influence the observable properties of galaxies. In particular, the dependence on "hidden variables" 
would be expected to increase the scatter in relations such as the mass-luminosity relation. The scatter seen here decreases 
strongly with mass. However, one might suspect that if it were possible to keep the same spatial resolution while working in a 
much larger simulation box, the effects of large-scale structure would lead to additional scatter at the high-mass end as well. 
Definite conclusions in this direction will unfortunately have to await further improvements in computational resources. 

Because the problem contains "hidden variables," there is no way to avoid uncertainties and scatter in the predictions 
based on averaged quantities. However, at least with regard to star formation, it is possible to avoid or at least reduce systematic 
errors resulting from not resolving small-scale densities by incorporating the available information into a characteristic star 
formation timescale t*. It can justifiably be argued that (for a given level of computer resources) some other hydrodynamical 
method such as SPH could achieve higher resolution in high-density regions than our methods. (I.e., some of the small-scale 
information which is "hidden" to us would thus no longer be hidden.) This better resolution has the potential to produce 
a more accurate estimate of cooling. Of course, at still smaller scales cooling is complicated by processes such as thermal 
instabilities, clouds, and the effects of a multicomponent medium. Nonetheless, it is still a valid question to ask how much 
influence fluctuations in density on scales below the resolution of our Models I - 10 could have. Although it is too early for 
a comprehensive answer, we can get a partial answer by comparing the results of the simulations with cooling enhancement 
parameter C = 1 (Models 3,7,9) with the C = 10 cases. The differences are not nearly as striking as the dependence on the 
supernova feedback parameter A. The reason may be that (enhanced) cooling is embedded in a feedback loop which couples 
all the processes and their timescales. If the timescale for an effect is already the shortest in a loop, then it makes little 
difference how short the timescale is. 

Moreover, a good characterization of the diffuse, lower-density, but hot gas is absolutely essential for a proper under- 
standing of the supernova feedback loop: First of all, some supernovae (about 3 %) are known to explode in the halos of 
spiral galaxies, so that the hot gas density is estimated correctly by our code. Second, even if a supernova explodes within 
the disk of a galaxy, the resulting supernova remnant will certainly deposit a significant proportion of its energy in the hot 
interstellar gas. Within the disk, which we do not resolve, the typical density and temperature are roughly 10~ 2 ' 5 cm -3 and 
10 5 ' 7 , respectively, and the cooling time would be on the order of 10 6 yr (McKee & Ostriker, 1977), which we underestimate 
somewhat even with enhanced cooling. However, a significant fraction of the supernova energy would escape to the halo, where 
the cooling rate is either correctly estimated or even somewhat overestimated by our prescription. Hence, it is impossible to 
argue that energy input due to supernovae is entirely locally dissipated and thus irrelevant for the gas dynamics on larger 
scales (Baron & White 1987; Navarro & White 1993; Kauffmann, White & Guiderdoni, 1993; Metzler & Evrard, 1994) 

In addition to the simulations reported in detail here, our experience with a large number of test simulations shows 
that the star formation history of galaxies and thus their observable properties such as colors and luminosities may vary 
quite strongly from one simulation to another depending crucially on the dynamics, especially interactions producing shocks. 
For example, if a galaxy merger just manages to occur at low redshift (shortly before the end of the simulation), then the 
ensuing burst of star formation dramatically increases the luminosities of the affected galaxies and shifts them toward "bluer" 
color indices. Hence, as a general statement concerning the interpretation of simulations we conclude that there is a very 
rich variety of interactions between large-scale and/or medium-scale structure and galaxy evolution that will require many 
numerical experiments and/or larger box sizes for an understanding. The trend toward increasing computer resources will 
allow larger boxes to be achieved and will thus hopefully lead to a more general understanding of statistical relationships 
between the environments of galaxies and their observable properties. 

Acknowledgements: 

We are very grateful to G. Bruzal for kindly making his models available to us. GY would like to thank C. S.Frenk, 
J.F.Navarro and A. Campos for very stimulating discussions. REK is grateful to V. Miiller, J. Miicket, and S. Gottlober for 
numerous productive discussions and moral support. GY also wishes to thank C. S. Frenk for inviting him to Durham. The 
visit was financed by the European Union Capital and Mobility Network Ulysses, whose support is gratefully acknowledged. 
This research was supported in part by NSF grant AST-9319970 (AK) and DGICyT of Spain (GY, project PB90-0182). REK 
was supported in part by a Fellowship (Kall81/1-1) from the DFG (Germany). 



REFERENCES 

Babul, A. & Rees, M.J., 1992, MNRAS, 255, 346 
Balbus, S.A. 1986, ApJ, 303, L79 

Bardeen, J.M., Bond, J.R., Kaiser, N., & Szalay, A.S., 1986 ApJ, 304, 15 
Baron, E., & White, S.D.M. 1987, MNRAS, 322, 585. 
Begelman, M.C., & McKee, C.F. 1990, ApJ, 358, 375 
Bertin, E. & Dennefeld, M., 1996, A&A, in press. 
Boris, J. P., & Book, D.L., 1973 J Comp. Phys., 11,38. 



Hydrodynamical simulations of galaxy formation: effects of supernova feedback 15 



Boris, J. P., & Book, D.L., 1976 J. Comp. Phys., 20, 397. 
Bruzual, G., 1983, ApJ, 273, 105. 
Bruzual, G & Chariot, S. 1993, ApJ, 405, 538. 
Buzzoni, A, 1989, ApJS 71, 817. 

Campos-Aguilar, A., Moles, M. & Masegosa, J., 1993, AJ, 106, 1784. 
Cen, R., 1992 ApJS 78, 341. 

Cen, R.Y. ,Jameson, A., Liu, F. & Ostriker, J. P., 1990 ApJ ,362,L41. 
Cen, R. & Ostriker, J. P. 1992, ApJ, 393, 22. 
Cen, R. & Ostriker, J. P. 1993, ApJ, 399, L113. 
Chariot, S. & Bruzual, G. , 1991, ApJ, 367, 126. 

Cole, S. Aragon-Salamanca, A., Frenk, C.S., Navarro, J. F. & Zepf, S. E., 1994, MNRAS, 271, 781. 
Colclla, P., & Woodward, P.R. 1984, J.Comp.Phys., 54, 174. 
Cowie, L.L., McKee, C.F., & Ostriker, J. P. 1981, ApJ, 247, 908. 
Dairies, S. J., Fabian, A.C., & Thomas, P.A., 1994, MNRAS, 268, 1060. 

de Vaucouleurs, 1977, in "The Evolution of Galaxies and Stellar Populations" eds. B.M Tinsley and R.B. Larson, Yale University 

Observatory. 
Dekel, A. & Silk, J., 1986, ApJ, 303,39. 
Ellis, R. S. et al; 1996, MNRAS, 280,235. 
Evrard, A.E., 1988 MNRAS, 235,911. 
Fall, S.M., Rees, M.J., 1985 ApJ, 298,18. 
Field, G.B., 1965 ApJ, 142, 531. 

Field, G.B., Goldsmith, D.W., & Habing, H.J., 1969 ApJ 155, L49. 
Field, G.B. & Saslaw, W.C.,1965, ApJ, 142,568. 

Gallcgo,J., Zamorano, J., Aragon-Salamanca, A., & Rego, M., 1995, ApJ, 455, LI 
Giroux, M. & Shapiro, P. 1986, ApJS, 102, 191. 

Hensler, G., & Burkert, A. 1990, Proc. " Windows on Galaxies," G. Fabbiano, et al. , Kluwer, Dordrecht, p. 321. 
Guiderdoni, B. & Rocca-Volmerange, B., 1987, A&A, 186, 1. 
Guiderdoni, B. & Rocca-Volmerange, B., 1990, A&A, 227, 362 
Hernquist, L., & Katz, N., 1989 ApJS ,70,419. 

Kang, H., Ostriker, J. P., Cen, R., Ryu, D., Hernquist, L., Evrard, A.E., Bryan, G.L., & Norman, M.L. 1994, ApJ, 430, 83 

Kates, R.E., Kotok, N., & Klypin, A.A., 1990 A&A 243,295. 

Katz, N. 1992 ApJ, 391, 502. 

Katz, N., & Gunn, J.E., 1991 ApJ, 377,365. 

Kauffmann, G., White, S.D.M., & Guiderdoni, B., 1993, MNRAS, 264, 201 
Klypin, A.A., & Kates, R.E., 1991, MNRAS, 251, 41p. 

Klypin, A. A., Kates, R.E., & Khokhlov, A., 1992 in 'New Insights into the Universe', eds. V. Martinez, M. Portilla, & D. Saez, Lecture 

Notes in Physics, Springer- Verlag, 171. 
Lada, C.J. 1985, in Star Forming Regions, eds. M. Peimbert & J. Jugaku, IAU Symp.115, Reidel, p.l. 
McKee, C.F. & Cowie, L., 1975, ApJ, 195, 715. 
McKee, C.F. & Cowie, L., 1977, ApJ, 215, 213. 
McKee, C.F. & Ostriker, J.P., 1977 ApJ, 218, 148. 
Metzler, C & Evrard, A. E., 1994, ApJ. 437, 564. 
Mihos, J.C. & Hernquist, L. , 1994a, ApJ, 437, 611. 
Mihos, C. & Hernquist, L., 1994b, ApJ, 431, L9. 
Miller, G.E., & Scalo, J.M., 1979 ApJS, 41,513. 
Miicket, J. P. & Kates, R., 1996, AIP preprint. 

Miicket, J. P., Petitjan, P., Kates, R, & Riediger, R., 1996, A&A , 308, 17. 
Navarro, J,F., & White, S.D.M., 1993, MNRAS, 265, 271 

Navarro, J,F., & Steinmetz, M.,1996, ApJ (submitted) astro-ph preprint 950643. 

Nulsen, P. E. J., 1986, MNRAS, 221, 377. 

Nulscn, P. E. J. & Fabian, A. C. 1995, MNRAS, 277, 561. 

Oran, E.S., & Boris, J. P. 1986, Numerical Simulation of Reactive Flow, Elsevier, New York 

Peebles, P.J.E. 1980, "The large-Scale Structure of the Universe" , Princeton University Press, Princeton 

Petitjean, P., Miicket, J, & Kates, R. ,1995, A&A, 295, L9. 

Raymond, J.C, Cox, D.P., & Smith, B.W., 1976 ApJ, 204, 290 

Shandarin, S.F., 1980 Astrophysica, 16, 439 

Sod, G.A., 1978 J. Comp. Phys., 21, 1 

Spitzer, L. 1990, ARA&A, 28, 71 

Tinsley, B. M., 1972, ApJ, 20, 283 

Thomas, P.A., 1988a, MNRAS, 235, 315. 

Thomas, P.A., 1988b, Proceedings of the NATO ASI "Cooling Flows and galaxy formation". Ed. A. C. Fabia, Kluwer, p235. 
Wolf, M.T., & Durisen, R.H. 1987, MNRAS, 224, 701 

Yepes, G., Kates, R., Klypin, A. & Khokhlov, A., 1995, Proceedings of the XVth Recontres des Moriond "Clustering in the Universe". 

Ed. S. Maurogordato, et al, in press. 
Yepes, G., Kates, R., Klypin, A. & Khokhlov, A., 1996. Proceedings of the UIPM-ECN Conference "Mapping, Measuring and 

Modelling the Universe". Ed. M. J. Pons, V. Martinez and P. Coles. PASP. pg 125. 



16 G. Yepes, R. Kates, A. Khokhlov and A. Klypin 



Table 1. Parameters of the simulations 



Parameters 








MODEL 


# 












1 


2 


3 


4 


5 * 


6 


7 


8* 


9* 


10 


Box Size (Mpc) 


5 


5 


5 


5 


5 


1.5 


1.5 


1.5 


1.5 


5 


Formal resolution at z = 1 (kpc) 


19 


19 


19 


19 


19 


5.8 


5.8 


5.8 


5.8 


9.5 


Mass resolution for dark matter (10 6 Mq) 


4 


4 


4 


4 


4 


0.1 


0.1 


0.1 


0.1 


4 


Feedback parameter A (eq.8) 


200 


200 


200 


200 


200 


200 


200 


100 


100 


200 


Cooling enhancement factor C (eq.3) 


10 


10 


1 


10 


10 


10 


1 


10 


1 


10 


Overdensity limit for star formation 


100 


100 


100 


50 


50 


100 


100 


100 


100 


100 


^luminous at 2 = (%) 


2.3 


2.4 


2.1 


2.7 


2.4 


2.1 










^luminous at Z = 1 (%) 


0.84 


0.77 


0.77 


1.4 


1.4 


0.68 


0.68 


0.27 


0.26 


0.78 


Z l/2 


0.65 


0.75 


0.63 


1.0 


1.2 


0.50 











Table 2. Volume average observational quantities at Z 
Quantity 

Luminosity density in B (10 7 Lq/Mpc 3 ) 
Bolometric M/L bol (Mq /Lq) 
SFR density (Mq /yr/Mpc 3 ) 
SFR density since Z = O.O5(M /yr/Mpc 3 ) 
Star density (10 9 M /Mpc 3 ) 





MODEL # 




1 


2 


3 


4 


5 


9.4 


8.0 


11.0 


8.5 


5.8 


31 


31 


24 


28 


85 


0.12 


0.13 


0.11 


0.15 


0.13 


0.08 


0.07 


0.06 


0.09 


0.08 


1.6 


1.6 


1.5 


1.9 


1.7 
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FIGURE CAPTIONS 

Fig 1 The density (open circles) and the pressure (filled circles) for a one- dimensional shock tube. The solid line represents the 
analytic Sod solution. 

Fig 2 The evolution of the density for a strong "2-D" cylindrical explosion iexplosion in an ideal gas with 7 = 5/3. A 3-D 
version of the code on a 100 x 100 x 1 grid was used with periodic boundary conditions along each spatial direction. The 
explosion energy is initially deposited inside a small cylinder of radius 2.5 cell widths located coaxially to the z-axis with 
P — 10 4 . In the left top panel, the density of every cell in the simulation is plotted against its distance from the center 
of the explosion. The top right panel compares the shock expansion with the analytical Sedov solution. The lower panels 
show the density at different moments. 

Fig 3 Effects of cooling and supernova feedback simulated in a single cell. No gravity and mass flow are taken into account 
(constant density). The time evolution of the mass fractions of cold gas, hot gas and stars are shown for different values 
of the supernova feedback parameter A: A = 1 (dash - dot), A = 20 (solid) and A = 100 (dashed). The gas, which has 
solar composition and initial temperature T — 10 6 K and density nu = 6 x 10~ 4 cm~ 3 cools with enhancement factor 
C = 10 (see text). 

Fig. 4 (a) A 3D view of the dark matter distribution in Model 1 at 2 = 0. A random sample of 10% of the dark matter particles 
is shown. Dimensions are in Mpc (Ho = 50 Km s _1 Mpc~ ) 

(b) Same as (a) but for the stars created in the simulation. The points in this figure do not directly represent the mass 
distribution of stars because star particles have variable masses. 
Fig. 5 The same as Fig. 4, but for Model 2. 

Fig. 6 (a) The distribution of temperature and densities in a slice of 1 cell size (39 kpc) passing through the center of the 
brightest galaxy of Model 1 at z = 0. Contours for quantities in parentheses are represented by thick solid curves. Lower 
contour values are represented by dashed curves. The density of the dark matter is expressed in units of the critical 
density p cr . The densities of the luminous matter ("stars") and the hot gas are expressed in units of the mean baryon 
density fltPa- 

(b) Same as (a), but in this case for Model 6 at z = 1. A cross section of the main filament running across the simulation 
box is clearly seen in the density contours of dark matter, gas and stars. The coordinate axes are scaled in proper 
(physical) units. 

(c) Same as (b) but for the high-resolution (256 3 grid cells) Model 10. 

Fig. 7 Relation between the luminous and dark-matter densities (both in units of the critical density). The density of the 
luminous matter and the density of the dark matter are plotted for every cell in the simulations, (a) Left panels are for 
Models 1,3, and 4 (bottom to top); Right panels are for Models 6,7, and 8 (top to bottom), (b) Same as (a) but for the 
high-resolution Model 10 at selected redshifts. 

Fig. 8 Fraction of luminous (star) mass to total mass as a function of the total mass for the galaxies in Models 1-4. Different 
symbols are used to plot galaxies at different relative distances to the most massive galaxy in the simulations. The dotted 
line represents the initial average baryon fraction (0.08). 

Fig. 9 Absolute visual magnitude of the galaxies in Models 1-4 versus total mass. Symbols as in Fig. 8. 

Fig. 10 UBV color-color diagram for the bright galaxies Mv < —15 in Models 1-4 (see figures). The long dashed curve represents 
the UBV position for main sequence stars of luminosity class V. The short dashed line is for galaxies (see text for details). 

Fig. 11 Star formation rates as a function of redshift for 3 different galaxies in Model 1 (top panel) and Model 6 (bottom panel). 
The solid line corresponds to the brightest galaxy in the respective simulation. 

Fig. 12 Dependence of star formation rates on physical parameters. In the upper panel, the upper curves are for the most massive 
galaxy, while the lower curves are for the least massive luminous galaxy. The solid lines correspond to Model 1 (also 
seen in Fig. 11), the dotted lines are for Model 3 (no enhanced cooling), while the dashed line is for Model 4 (lower 
star-formation threshold). The lower panel is for high-resolution simulations: the solid lines are for Model 6, the dotted 
lines for Model 7 (no enhanced cooling), and the dashed lines for Model 8 (lower supernova feedback parameter). 



